#include "EBFV1_elliptic.h"
namespace PRS
{
	int EBFV1_elliptic::divergence_G(Mat G, int id0, int id1, double MobIJ, double aux){
		double Gij_I[2] = { MobIJ*aux,-MobIJ*aux};
		double Gij_J[2] = {-MobIJ*aux, MobIJ*aux};
		int rows[2] = {id0-1,id1-1};	// where Gij must be assembled in global G
		int cols[2] = {id0-1,id1-1};	// where Gij must be assembled in global G
		ierr = MatSetValues(G,1,&rows[0],2,cols,Gij_I,ADD_VALUES);CHKERRQ(ierr);
		ierr = MatSetValues(G,1,&rows[1],2,cols,Gij_J,ADD_VALUES);CHKERRQ(ierr);
	}

	int EBFV1_elliptic::divergence_E(Mat E, int dim, int id0, int id1, double MobIJ, EdgeData* edata){
		int i,j;
		double Eij_part1[3] = {.0,.0,.0};
		for (i=0; i<dim; i++){
			for (j=0; j<dim; j++){
				Eij_part1[i] += edata->Cij[j]*edata->aux_Matrix_E[dim*j+i];
			}
			// if an edge is on partition boundary, the result from all processors contribution for this edge must be equal as if was ran in serial
			// Eij (serial) = Eij (parallel) if edge is on partition boundary
			Eij_part1[i] *= MobIJ;///(edata->numRC+1.0);
		}

		// fill Eij matrix
		double Eij_row1[2*dim], Eij_row2[2*dim];
		for (i=0; i<dim; i++){
			Eij_row1[i] = Eij_part1[i];
			Eij_row1[i+dim] = Eij_row1[i];
			Eij_row2[i] = -Eij_row1[i];
			Eij_row2[i+dim] = -Eij_row1[i];
		}

		// where Eij must be assembled in global E matrix
		int pos1 = dim*(id0-1);
		int pos2 = dim*(id1-1);
		int idxn[2*dim];
		int idxm[2] = {id0-1,id1-1};
		for (i=0; i<dim; i++){
			idxn[i] = pos1+i;
			idxn[dim+i] = pos2+i;
		}
		ierr = MatSetValues(E,1,&idxm[0],2*dim,idxn,Eij_row1,ADD_VALUES); CHKERRQ(ierr);
		ierr = MatSetValues(E,1,&idxm[1],2*dim,idxn,Eij_row2,ADD_VALUES); CHKERRQ(ierr);
		return 0;
	}

	int EBFV1_elliptic::gradient_F_omega(Mat F, int dim, int id0, int id1, EdgeData *edata){
		VertexData* v_0 = edata->vertex_0;
		VertexData* v_1 = edata->vertex_1;
		double nrc = 1.0;//(double)edata->numRC + 1.0;
		double volumeI = v_0->volume;
		double volumeJ = v_1->volume;
		double aux1 = (1./(2.*volumeI))/nrc;
		double aux2 = (-1./(2.*volumeJ))/nrc;

		// Assembly Fij
		double Fij_column1[2*dim], Fij_column2[2*dim];
		for (int i=0; i<dim; i++){
			Fij_column1[i] =  edata->Cij[i]*aux1;
			Fij_column2[i] = Fij_column1[i];
			Fij_column1[i+dim] =  edata->Cij[i]*aux2;
			Fij_column2[i+dim] = Fij_column1[i+dim];
		}
		int pos1 = dim*(id0-1);
		int pos2 = dim*(id1-1);
		// says on which rows Fij must be assembled into Fg
		int idxm[2*dim];
		for (int i=0; i<dim; i++){
			idxm[i] = pos1+i;
			idxm[dim+i] = pos2+i;
		}
		// says on which columns Fij must be assembled into Fg
		int idxn[2] = {id0-1,id1-1};

		ierr = MatSetValues(F,2*dim,idxm,1,&idxn[0],Fij_column1,ADD_VALUES);CHKERRQ(ierr);
		ierr = MatSetValues(F,2*dim,idxm,1,&idxn[1],Fij_column2,ADD_VALUES);CHKERRQ(ierr);
		return 0;
	}
}
